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Abstract 

Many models have been developed to study the role of branching 
actin networks in motility. One important component of those models 
is the distribution of filament orientations relative to the cell mem- 
brane. Two mean-field models previously proposed are generalized 
and analyzed. In particular, we find that both models uniquely select 
for a dominant orientation pattern. In the linear case, the pattern 
is the eigenfunction associated with the principal eigenvalue. In the 
nonlinear case, we show there exists a unique equilibrium and that the 
equilibrium is locally stable. Approximate techniques are then used to 
provide evidence for global stability. 

1 Introduction 

Actin is a protein involved in many cellular processes ranging from regulating 
gene transcription to acting as a motor in cell motility j?]. It is one of the 
most conserved proteins and is present in almost all eukaryotic cells. Actin 
monomers polymerize into thin filaments which form highly branched net- 
works near the leading edge of motile cells [16j. While actin monomers will 
spontaneously polymerize in physiological conditions, inside these branched 
networks, new filaments are generated by branching off of existing filaments 
[19]. New filaments are nucleated by the actin related proteins 2 and 3 com- 
plex (Arp2/3). To maintain a consistent supply of actin monomers, actin 
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filaments are eventually severed and depolymerized. Filament density is reg- 
ulated by capping protein binding to the filament tips, ceasing polymeriza- 
tion [5j. Combined with filaments growing by the addition of new monomers, 
these processes create a dynamic network that serves as the engine in certain 
types of cell motility UHl ED- 

Any individual filament in an actin network can be partially character- 
ized by the angle between it and the normal direction of the membrane. One 
obvious question is whether or not these angles form any regular pattern. 
While the question has not been extensively studied experimentally, there is 
some evidence that the networks indeed organize into regular patterns rela- 
tive to the cell membrane [T^ |221 ESI EZ]- A few models have been proposed 
to explain the existence of such patterns [T^ EH EH]- While these models 
have been numerically studied, there has been no rigorous work proving the 
existence, uniqueness or stability of these solutions. This article presents 
a few results that characterize the solutions to two equations modeling the 
angular density of branching actin networks. 

All of the models proposed to explain the orientation distribution have 
used a continuum approximation. There is some question as to whether or 
not ignoring the stochastic fluctuations of actin networks is justified |21] . 
However, none of the models make specific predictions about the kinetics 
of network organization, and there is currently no evidence that correlations 
between filaments lead to changes in the equilibrium orientation pattern. For 
the rest of this article, we will assume the approximation is justified and focus 
on long-time equilibrium behavior. 

Some of first few models to study orientation patterns in actin and simi- 
lar networks studied the existence and persistence of peaks in the orientation 
pattern |9l [131 E] • The analysis was based on Fourier series and small per- 
turbations which greatly limited their generality. Their analysis led to the 
qualitative result that peaked orientation patterns are likely to be observed. 
Similar methods have been used on models of orientation and space [U E]. 
Stability analysis has also been done on similar models, termed "ring models" 
in the neuroscience literature [3] [28]. 

The first model we consider here was proposed by Maly and Borisy |12] . 
Their insight was that if filaments were capped at different rates based on the 
filament orientation, filament branching and capping could generate stable 
orientation patterns. The model they proposed takes into account branching 
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and capping explicitly and filament growth implicitly. New filaments branch 
off of existing filaments at a characteristic angle ~ 70° with some variance 
around that. We can write out a branching kernel as a probability of a mother 
filament with angle 6m having a branched daughter filament with angle Oe,: 

i3(0)=P(^^D = ^A/-0) (1.1) 

Adding up the contribution of all filaments with density gives the total 
branching rate at angle 9: 

BK{0) (X j B{d - (t))u{(l))A(t) (1.2) 

They also proposed that the capping rate was proportional to the amount of 
time the filament would be not in contact with the leading edge, much like 
[15], but used the capping function ^^J^g^ . The Maly and Borisy model only 
considered filaments growing faster than the leading edge, i.e. filaments with 

orientation 1^1 < ^crit = arccos where is the velocity of the leading 

edge relative to the maximum velocity of filament growth. Combining the 
two terms gives the full equation: 

u{9,t) = \ I B{9 - 4>H4>) - (1.3) 

where ii indicates the time derivative. The equation is defined on (— ^crit, ^crit) x 
]R+ with absorbing boundary conditions. 



Maly and Borisy performed two analyses on (1.3). The first analysis was 



to approixmate solutions of (1.3) by solving the equation for two points in 
orientation space. Solutions to the two-point approximation supported the 
argument that the equation selected for a unique orientation 'type' that grew 
exponentially at the fastest rate. The second analysis was to use numerical 
quadrature [2j to approximate the eigenfunctions of the right-hand side of 



(1.3). However, the existence and uniqueness of the eigenfunction solutions 



were never rigorously shown. They explained their results by using an evo- 



lutionary selection metaphor. In this article, we show that a version of (1.3) 
with stricter hypotheses on the capping rate uniquely selects for a most 'fit' 
orientation pattern with a fitness function defined on the unit ball of orien- 
tation functions. 
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A very similar model was proposed by Weichsel and Schwartz to 
explain both the orientation patterns and the velocity of a growing actin 
network pushing against a given force. There were two primary differences 
between their model and the Maly and Borisy model. First, orientations 
were defined on the entire circle, 5*^. The second difference was to normalize 
the total branching rate to the constant A, which ensures that solutions have 
bounded total density. The Weichsel and Schwartz model is: 

jg,u{(p,t) d(f) J 

where the capping rate is a constant plus a term proportional to the differ- 
ence between the velocity of the leading edge and a filament with a given 
orientation: 

K{e) = k + C{VLE-Vocos{e))+ (1.5) 

where Vle is the velocity of the leading edge, Vq is the rate of filament growth, 
and {x)~^ is the positive part of x. 

Weichsel and Schwartz performed the same two analyses as in the Maly 
and Borisy paper. They found that, for certain parameters, there were two 



equilibria in the two-point approximation to (1.4), where one equilibrium is 
stable and the other is a saddle. Finally, they used numerical techniques 
to calculate the equilibrium distributions. The results in this article explic- 
itly contradict their assertion of multiple equilibria, but they show the local 
stability of a unique, positive equilibrium. However, it is important to note 
that the equilibrium is unique for a fixed B and k. Changing the load force, 
concentrations of branching and capping proteins, or other experimental ma- 
nipulations could change the structure of the unique equilibrium. 



1.1 Overview of Results 

The two equations we specifically analyze are: 

u{e, t) = x(^B^uj {6, t) - k{6)u{6, t) (1.6) 

u{e, t) = \ (b ^ u) {9, t) - k{9)u{9, t) (1.7) 
J u[ui, t)aui \ / 
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Equation 1.6 is our generalization of the Maly and Borisy [12] model, and 
equation 1.7 is our generalization of the Weichsel and Schwarz [26] model. 
Both equations are defined on the circle with B > being the branch- 
ing kernel which generates new filaments and k{6) > being the variable 
capping rate which eliminates filaments. The hypotheses on each function 
are relatively weak. Without loss of generality, we will assume that B has 
unit integral (/ B{u) du = 1). We will also assume that B is continuous and 
even. Finally, we will assume that k is a-Holder with a = 2. There are other 
combinations of hypotheses that lead to similar results. A paper studying 
necessary and / or sufficient hypotheses for the existence of eigenvalues for a 
generalized nonlocal diffusion equation is forthcoming. 

In agreement with the paper [12], the first-order branching equation 
uniquely selects for an optimal orientation pattern. However, the Weich- 
sel and Schwarz paper [26] suggests that there might be multiple equilibria. 
We show that the zeroth-order branching equation also uniquely selects for 
a unique equilibrium orientation pattern. 

The first two results characterize solutions to the first-order branching 
equation (1.6). Theorem [T] shows that the spectrum of the operator defining 
(1.6) is dominated by an isolated, simple principal eigenvalue with strictly 
positive eigenfunction. While that eigenvalue may be positive or negative in 
general, long-time solutions to (1.6) are therefore dominated by the exponen- 
tial increase or decay of the principal eigenfunction. Proposition [3] says that 
for given B and k, there exists only one A such that (1.6) has a non-trivial 
equilibrium. 



The rest of the article is dedicated to analyzing (1.7). Proposition 0] gives 



the existence and uniqueness of a non-trivial equilibrium. Linear stability 
analysis combined with Theorem |6] implies that the equilibrium is locally 
stable. Finally, numerical simulations and a perturbation analysis are per- 
formed to provide evidence that (1.7) is globally stable. 



2 First-order Branching 



The first result uniquely characterizes the dynamics of (1.6). Define A to be 



the linear operator on the right-hand side of (1.6). It is easy to observe that 



A is resolvent positive and therefore generates a positive, Cq semigroup. The 
long-time dynamics are described by Theorem [T} 
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Theorem 1. A has an isolated, algebraically simple principle eigenvalue with 
positive eigenf unction. 



Proof. It is straightforward to observe that A is a bounded hnear operator 
on both of the spaces L}{S^) and L'^{S^). In fact, A is self-adjoint on L^. 
The convolution B-kui?, also compact. To see that, take a sequence {"UjjjeN 
with IImjIIi = 1. We first have that the sequence B-kUj is uniformly bounded: 



\B -k UjWoo < \\B\ 



= Moo (2.1) 

by Young's inequality. We can also show that the sequence is equicontinuous: 



sup \{B -k Uj){x) 
\x-y\<S 



iB^u,){y)\ 



sup 

\x-y\<5 



{B{x - z) - B{y ~ z))uj{z) dz 



< sup \B{x) 

\x-y\<5 



Uj{z) dz 



(2.2) 



where h{5) — )■ as 5 — )■ since B is uniformly continuous. Arzela-Ascoli 
ensures the existence of a convergent subsequence and a continuous limit 
function. 

Next, we will characterize the spectrum of A outside of the negative range 
of K : —R{k,) = [— sup K, — inf k]. Observe that we have the two equivalence: 



B-k V — KV - 

B^v 



liv = f 



f 



K, + fi K + /i 

Since K + /i is bounded away from zero. 



(2.3) 



-f- G whenever f & L^. ^ 
is a compact operator, so by the spectral theory of compact operators, (2.3) 
is solvable unless 1 is an eigenvalue of Given the equivalence, we can 



conclude that the residual spectrum of A, if it is non-empty, is contained in 
—R{k,). Now, assume fi is in the continuous spectrum of A. That implies 
there exists {uj}, \\uj\\i = 1, where: 



lim \\Auj — /iMjIli = 



By the compactness of B-k, there exists uj^, with B-kUj^ — >G C{S^): 



= lim WB-ku 



Ik 



lim \\w 

k—>oo 



l^)Uj, 



jk II 1 



(2.4) 



(2.5) 
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The above limit implies that -Uj^. — )■ and that G C{S^) is and eigen- 
function of A with eigenvalue /x. The two preceeding results imply that all 
elements of the spectrum outside of —R{k) are eigenvalues. In fact, we have 
also shown that the spectrum of A on is equal to the spectrum on out- 
side of —R{k). Thus, we can construct the remainder of the proof defining 
A on the Hilbert space with the usual inner product. 

It remains to show that there exists an element of the spectrum /iq > 
— inf K. The technique presented in |1U| would generate this conclusion under 
the hypothesis that B{0) > 0. However, that is not physically realistic as new 
filaments will always branch off of an existing filament with an orientation 
different that the mother filament. In order to drop the hypothesis that 
B{0) = 0, we add the hypothesis that n is a-Holder continuous with a = 2. 
Since A is self-adjoint, we know the largest element of the spectrum over 
is: 

yUo = sup {Au,u) (2.6) 
||m||2=i 

A is also compact over L^, so /io is an eigenvalue as long as /iq > — inf k. It 
suffices to show there exists u & L? such that: 

(Am,u, ) + inf > (2.7) 

By the continuity of n and compactness of the circle, k,(9) — inf k = for at 
least one 9. For the sake of notation, define g{9) = k,{6) — inf k. 

Observe that for any function of the form u{6) = c + f{6) > with 
c, f{e)>0 where J^, u{e) = 1: 

{Bi.u,u) = jj B{B -uj){c+ f{uj)) dw(c + f{e)) AO 
>c[[ B{9 -uj){c + f{e))dLod9 = c 



since f{6) > and J B = 1. Without loss of generality, assume g{0) = 0. By 
our hypothesis that k is Holder continuous, we have the inequality: 

9(0) = 9(9) - giO) < Q\e - 0\' = Q9' (2.8) 

Define as: 

_ l[-s,s] 

2e 
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where 1 is the usual indicator function. Note that J f — 1. We now have 
the two relations: 



51 -£ 



2s^ 



£ e 

fe{efg{e) de = J fe{efg{e)de <^jQe'de = g| 

51 -e 

Combining those relations gives: 

(5(c + A), (c + /,)) = J g{e){c + uo)yde 

51 

= J g{9) (c^ + 2cMe) + Udf) de < Rc' + cQ^ + Qj 
51 



where R^J g(e) dO. 

First, assume i? > 27r and Q >1. Fix c = jji, £ = and c' — 1 — 27rc < 
1. Fix = c + c'/e- Putting all of the above together gives: 

{Au + (inf k)u, u) = {Bi.{c + c'f,),c + c'/s) - {9(0 + c'/.), c + c'/.) 
> c - i?c^ - cc'g^ - c"Qj 

- 2R 4i?2 R s ^3 
^ 1 Q 1 Q 



AR 3R SR^Q'^ 6RQ 
1 1 11 
- 4R ~ 24R ~ 6^ ~ 24R ^ ° 

If i? < 271, then set c = 2^ and c' = 0. If i? > 27r and Q < 1, set c and 
c' as above and £ = That shows we have constructed such a m and 
max (72 (A) > — infre. 

Now that we have a principal eigenvalue, to show that it is isolated, take 
a sequence /ij — > /Iq with associated eigenfunctions Uj. Choose a subsequence 
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such that B-kUj^ is convergent: 



= hm \\B-kUj^ II2 

fc->oo 

= hm II (k + /ijjuj, - (k + II 2 

fc— s>oo 

= hm ||(/t + /io)Kfe -«ife+i)ll2 

fe— >oo 

< hm (inf k + iiQ)\\uj^ - Uj \\2 = (inf k + no)V2 > (2.9) 

fe— 5>00 

since Uj and u^,, are orthogonal and fiQ > — inf k. 

The last thing to show is that the eigenvalue is simple with positive 
eigenfunction. Recalling the equations in (2.3), we know that the following 
eigenvalue equation has at least one solution: 

u (2.10) 



A'u = -^fjf^ is obviously a positive and compact operator on the Banach 



B*u 

K + flQ 

space of continuous functions. The Krein-Rutman theorem implies that A' 
has an eigenvalue equal to its spectral radius. Assume that spectral radius 
P(A') > 1: 

—— = P(A')« ^ 

K + /io 

i3 tIt M — (km + /io)M = (p(A') — 1)(k + /io)^ <^=^ 

Au - fiou = {p{A') - 1){k + ij,o)u (2.11) 

That last equality implies: 

(Am - /ioM, u) = (p(A') - 1){{k + /io)u, u) > (2.12) 

in contradiction to our definition of fiQ. Therefore, p(A') = 1, and there 
exists a positive u that solves (2.10). 

The last thing that remains to be shown is that the eigenvalue is simple. 
The Krein-Rutman theorem implies that it is sufficient to show that C^'^ 
is strongly positive for some n. That is shown Lemma [2j With Lemma ^ we 
have shown that (^^'^ is a strongly positive operator with leading eigenvalue 

1. Since A' and (^A'j have the same eigenvalues, (^A'j having a simple 
leading eigenvalue implies the leading eigenvalue of A' is also simple. □ 



9 



Lemma 2. For any given kernel B, there exists n such that ^A' j is strongly 
positive, i.e.: 

u>0 ^ [^'Yu > 
whenever u> is a continuous function not uniformly zero. 

Proof. We are only concerned with whether or not ^A'j is positive and 

not on the specific value of ^A' j , so it is sufficient to show the result for 

B" where Bu — B -k u. Define E to be the cr— algebra associated with the 
Lebesgue measure on S^. We can define the set mapping T : E — > E as: 

TO = supp{-B^lf^} (2.13) 

where f2 G S. Choose some open interval {y — 6,y + 6) C supp{i3} for y ^ 0, 
y e M\Q. Assume Q contains an open interval {x — e,x + e). Observe that 
for every z & {x + y — S, x + y + S): 

x+e y+e 

Bi.la{z) = J B{z-s)ln{s)ds> J B{z-s)ds= J B(^{z-x+y)-s'^ ds' > 

x—e y—e 

(2.14) 

The above argument also holds for z e {x — y — S,x — y + S). Notice that 
while the existence of the interval was used in the above calculation, there is 
no explicit dependence on e beyond that s > 0. Iterating T, we can observe 
that: 

T^^O C Uo<,<n ((x + 2jy -5,x + 2jy + 5)U{x- 2jy -5,x- 2jy + 5)) (2.15) 

By the fact that x + 2y is an irrational rotation, {x + 2jy}j^^ is dense 
in and the sets {{x + 2jy — S,x + 2jy + S)}j^^ form an open covering 
of S^. The compactness of implies that there exists n E N such that 
C Ui<j<„(a; + 2jy — 6,x + 2jy + 6). Rotational symmetry in implies 
that n has no dependence on x. Fix some function u G C{S^) with m > 
and u not uniformly zero. We can set Q = supp{m} and observe ft contains 
an open interval containing some x'. The above discussion implies: 

T^'^n D Ui<j<n{x' + 2jy -6,x' + 2jy + 6) D (2.16) 

The above set relation implies B^"k > 0. □ 
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Another small proposition to characterize solutions to (1.6): 



Proposition 3. Given B and k, there exists precisely one A such that (1.6) 
has a stable, non-trivial equilibrium. 

Define the operator Ax as: 

A\U = XB -kU — KU 

Define the related operator and inner product spaces: 



Ki B-ku . „ , 



K 



We can now prove the result. 

Proof. Observe A' is self-adjoint with respect to (■,■)«• Also, we have the 
relation between A^ and A': 



(Aau, u) = jj ^^{^ - u{e) - Ku{ef dO 

B{e-uj 



x^—-^u(u)duu(9)K(e)d9- / u{eyK(e)de 

= X{A'u,u)^ - {u,u)^. 

From the proof of Theorem [T} we have that A' has a simple principal eigen- 
value, /i'o > 0. We know that /io and /Xq can be defined by the following: 

(A'm,m)«, 



(Am, u) 
/io = sup — — 



and /iQ = sup 

There is a relationship between the sign of /io and /ig: 

X{A'u,u, )«; - {u,u, )k 



(2.17) 



sign[A/io - 1] = sign 



sup 



sup 
sup 

Il«ll27^0 



Sign 



A(A'm,m, )k - {u,u,)^ 



sign 



sign 



{Axu,u) 
{Axu,u) 



sup 

||«||2^o {u,u) 



sign[/io] 



(2.18) 
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The argument holds since 



{u,u)k 

{u,u) 



> by hypothesis of the supremum and 



does not change the sign of the argument. We know that (1.6) has a non- 
trivial equilibrium if and only if Aa has a zero eigenvalue. That equilibrium 
is stable if and only if all of the other elements of the spectrum have negative 
real part. That is the case if and only if the zero eigenvalue is the largest 
eigenvalue, i.e. /io = 0. The calculation above therefore implies there exists 
only one A where (1.6) has a stable equilibrium since there is only on A such 
that fiQ = A/iQ — 1 = 0. □ 



3 Zeroth-order Branching 



Proving the existence of solutions to (1.7) does not require any sophisticated 



machinery. Define G{u) to be the nonlinear operator that defines the dynam- 



ics of (1.7). It is straighforward to show that G is locally-Lipschitz. Global 
existence for alH > can be shown by observing that solutions with positive, 
integrable initial data stay in the closed set: 

U := {u e : u{9) > and < c < ||m||i < c' < oo} 

G{u) is uniformly Lipschitz on U. It is obvious that solutions with positive 
initial data remain positive. Simply observe for any angle 6* with u{6*,t) = 
where u{-, t) > 0: 

u{e,t) = {Bi.u){e,t) - K{e)u{e,t) = > o (3.1) 

It is also straightforward to show that there exists c and c' for the definition 
of 14. First, observe that 

^ J u{e, t)de = Xo- J K{e)u{e, t) de (3.2) 

The mean value theorem gives: 

(inf/t) j u{9,t)d9< J K{0)u{0,t) d0 < {sup k) j u{9,t)d9 (3.3) 
We can then write out explicit expressions for c and c': 

c = min|y"n(^,0)d^, c' = max ^(0, 0) d^, (3.4) 
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All that remains is to show that G is uniformly Lipschitz on U. The derivative 
of G is equal to: 

DGM = - <BHO) - /.Md. /'^7";<;>, (3,5) 

It is easy to see that a coarse estimate for C is: 

2 

-DG^ op ^ ~1~ oo 
c 

which implies that G is Lipschitz on U. The standard Picard-Lindelof argu- 
ment is sufficient to show existence on x M"*". 

The first result is an existence result: 



Proposition 4. A function u E is an equilibrium of equation (1.7) if and 
only if it is a solution to the eigenvalue problem: 

^=,u(e) (3.6) 

where yU 7^ and J^^ u{u) du ^ 0. 

Proposition 111 is in contradiction to the hypothesis in [26] that there are 



multiple equilibrium solutions to (1.7). 

Proof. Assume you have an equilibrium u E with J u{uj) dw 7^ 0, i.e. 
G(m) = 0. We know that: 



/^i u{ijj) do; 
Simple algebra gives: 

{Bi.u){e) _ 



u{uj) du u{6) 



That implies u is an eigenfunction with eigenvalue J u{u) do; 7^ 0. For the 
other direction, assume that we have: 

= fiU 

13 



and the listed hypotheses above. Simple algebra again: 



B -kU — fJiKU = 

Assume J^^ u{uj) du = 1. This is justified as long as J^^ u{u) du ^ 0. 

The existence of at least one positive eigenfunction with non-zero integral 
is ensured by the Krein-Rutman theorem as in Theorem [T| f = ^qUq is now 



an equilibrium to (1.7). Lemma 2 ensures that m > 0. By the self-adjointness 



of A', we have the all eigenfunctions Uk 7^ Mq are orthogonal to Uq: 
= {uo,Uk)n = [ uoie)uki9)Ki9)de 



51 



However, we know that kuq > 0. That implies the above can only be true 
if = almost everywhere or Uk is negative on some set with non-zero 
measure. □ 



The next result implies stability of (1.7) in a local sense. Before the proof. 



a quick, basic lemma from complex analysis: 

Lemma 5. Assume x G M and < x < y. There exists a real function 
P{a, y) > for « 7^ G M such that: 

^ < 3^ (3-7) 



|x + za;| X + (3 

Proof. Direct calculation shows: 

1 1 

< 



|x + ia| X + (3 
{x + P'f <x^ + <^ 
< + - (x + 

0<a^- (3^- 2x13 (3.8) 

Choosing f3{a,y) = min |l, completes the proof. □ 

Theorem 6. The spectral bound of the linearization around uq, the equilib- 



rium of (1.7), is strictly less than zero. 
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Define the spectral bound to be: 

s(A) = sup{3?/i : fi G (t(A)} (3.9) 

We will show that the right half of the complex plane is contained in the 
resolvent. For these purposes, we will only consider A as an operator over 
L^, but the results are immediately generalizable to using the techniques 
in the proof of Theorem [T| For all 7 such that 9^7 > and 7 7^ 0, the 
proof will consist of constructing a Neumann-type series and showing the 
series converges to a bounded operator. The proof is completed by showing 
^ a{DGuo). 



Proof. First, fix 7 > with 7 G M. We can solve explicitly for the resolvent. 
Fix / in and assume there exists v such that: 



Av — jv = f 



(3.10) 



We will derive a Neumann-type series to show the existence of such a v. 
Expanding out A and rearranging gives: 



f 



K + 7 K + 7 
From previous results, it is easy to see that: 



(3.11) 



K + 7 



< r 



Bi^ 



where r(-) is the spectral radius. Define : 
series gives us the explicit form for (A — 7)^^: 



K+7 ■ 



j=0 



f 



K + 7 



(3.12) 

The usual Neumann 
(3.13) 



The convergence of the above series implies that (A — 7) ^ is a bounded 
operator and 7 G p(A), the resolvent. 



Equation (3.11) is equally valid for complex 7. Assume 3fJ7 > and 
7 7^ 0. In order to show that (3.13) still holds, we need to show that the 
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spectral radius is strictly less than one. Using Lemma [5} we can provide a 
bound on the numerical radius n(-): 



n{B^) = sup |(Bt>,t>)|= sup 



f 2 = 1 



\M2=1 

< sup 

||d!|2=i 

< sup 

||l'||2 = l 



\{B^vmv{e)\ 
\{i3i<v){e)v 



v{9) d9 

de 

de < 1 



(3.14) 



The Neumann series in equation (3.13) therefore converges and 7 G p(A) for 
all 7 7^ where 3^7 > 0. 

The last remaining case to prove is that G p{DGuo)- Assume G 
a{DGuo)- That would imply that: 



v{uj) did kuq 



(3.15) 



The nullspace of the left hand side is spanned by Uq. Since v = Uq does not 
solve the equation, we must have J f 7^ 0. However, the Fredholm alternative 
states that the above is only solvable if the right hand side is perpendicular 
to uq. As huq is not perpendicular to uq, the equation is not solvable. □ 



4 Approximate Methods 

The stability result in Theorem [6] is a local result, so approximate methods 



were used to characterize the global behavior of (1.7). First, a set of nu- 
merical simulations were performed that showed global exponential stability. 
Then, we present a perturbation expansion that provides evidence towards 
global stability along with a discussion of the limitations of the perturbation. 
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4.1 Numerical Simulations 



The two branching kernels were: 



and 

^ + if|^|<i 

82(9) oc < 4.1 

^ \0 if 1^1 > 1 

where B2 was normahzed to have integral one and 0(^) was a von Mises 
distribution: 

. 5^El^l£2!M (4.2) 

27r/o(cr ^) 

where o" = -Bi was based on the branching kernels used in [T2] and |26] . 
Those two papers used truncated Gaussian distributions centered around 
±70°. Here, the von Mises distribution was used to avoid truncating the 
Gaussian or using the more complicated, formally correct wrapped Gaussian 
distribution. Also, the offset of ±| was used to simplify the radians conver- 
sion. Finally, the constant a was chosen to be in line with previous numerical 
studies |23t [2^ 126] . Both branching kernels were and symmetric. 

The two capping functions used were: 

Ki{d) = l-]^ cos(^) and K2{9) = 1 + ^ cos(4^2) (4.3) 

The first capping function, Ki was based upon |26], and the second was chosen 
to have multiple minima and maxima and non-uniform oscillations. 

Finally, the four initial conditions chosen were: 

ui{e) = 1 U2{e) = - — ^ 



-/2 



3 I 



«3W = lr-,.) + l(_^,_^) u,{9) = ^^l^ (4.4) 

They were chosen to include a mix of symmetric, non-symmetric, smooth 
and non-smooth functions. Also, U2 was chosen so that U2 G (L^\L^). 



17 




Angle Angle 

Figure 1: The plots here show the equihbrium distributions calculated using 
the method in section [B| for the four systems. A) Bi and ki, B) Bi and k,2, 
C) B2 and Ki, and D) B2 and K2. 
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The first calculations run were to estimate the equilibrium distribution 
using the method in section |B} The method was iterated until the || ■ ||i differ- 
ence between successive iterations was less than double precision. The equi- 
librium distribution appeared to have a qualitatively stronger dependence on 
K than on B as can be seen in Figure [T| 



A) 



B) 




C) 



40 60 

Time 




D) 



40 60 

Time 




40 60 

Time 




40 60 

Time 



Figure 2: The observed L} distances from equilibrium as a function of time 
are plotted here on a log scale. The four lines on each plot each correspond 
to the solution starting with a different initial condition. A) Bi and Ki. B) 
Bi and K2. C) B2 and Ki. D) B2 and K2. 



All of the simulations run converged (asymptotically) exponentially to the 
equilibrium. The equilibrium was calculated by iterating A' as outlined in 
the following section. Figured shows the distance between the simulation 
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result and the calculated equilibrium on a log scale. The log scale was used to 
make the graphs legible and to show the exponential convergence. Moreover, 
all of the initial conditions appear to asymptotically converge at the same 
rate. That gives evidence that solutions are globally exponentially stable. 
The simulations provide evidence for exponential stability with some constant 
determined solely by B and n. 



4.2 Perturbation Expansion 

The last approximate technique we will use is to study the perturbation 
expanion of the zeroth-order branching equation. For this section, we will 
again consider functions u ^ L}. However, we will again appeal to Hilbert 
space techniques when necessary. 

Assume that the capping function can be written out as: 

n{e)=c + e(t){e) (4.5) 

where is smooth, has integral zero, and reasonably small so that £0 is close 
to zero. The equation of motion is thus: 

^) = V I - + (4-6) 

J u\uj, t) au 

For the Cauchy problem with u{9,Qi) = u*{6), we will look for solutions of 
the form: 

oo 

u{9,t) = Y,e'u,{e,t) (4.7) 

j=0 

with the initial conditions U(){6, 0) = u* and Uj{6, 0) = for all j > 1. Show- 
ing that both Mo and Ui converge to the equilibribium defined in equation 
(3.6) provides some evidence for the stability of the full equation. 



First, we will calculate the first couple of terms of the equilibrium using 



equation (3.6) 
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where we have the two expansions: 



w{e) = J2e'w,{e) /. = 5^5V, (4.9) 

j=0 j=0 



Expanding out the terms in (4.8) gives: 

d=o / \j=o / \i=0 / \i=0 

oo j (_A\k oo j 

E E ^ ^^-^ = E E (4-10) 

j=0 k=0 j=0 k=0 

The first term (e^) is simply the eigenproblem for the unperturbed problem: 

B-kWo = CfioWo (4-11) 

The only positive solution of the above is /io = ^ and wq constant. Since we 
know the equilibrium has the same integral as the eigenvalue, we know that 
Wo = The term gives: 



c (? 



l^oWi+jJiWo (4.12) 



Filling in the known quantities and rearranging gives: 

B^wi-wi = ^ + !^ (4.13) 
c 2% 

The left hand side is a self-adjoint Fredholm operator with nuUspace spanned 
by the constant function, Wq. By the Fredholm alternative, we know that 
(4.13) is solvable if and only if the right hand side is orthogonal to Wq, 
i.e. (j){6) + ^ has integral zero. To have integral zero, we know that fii = 
-^/0(a;)dI; = O. 

We can write out Wi in terms of the Neumann series. By hypothesis, (j) is 
bounded, and therefore (j) ^ L'^. We know that B-k has operator norm strictly 
less than one on the space of functions orthogonal to wq: {wq}^ C L^, which 
implies the Neumann series converges in norm: 

oo 

Wi = ^B^0 (4.14) 

3=0 

21 



where Bu = B -k u and W^^u = B -k Wu. The above expansion imphes that 
J Wi = 0. We know that j Bn = j u which gives: 

/n 
J2{^'(l>)id)d0 = (4.15) 
j=0 

Since norm convergence imphes weak convergence, we can finish the proof 
by observing: 

hm f wi- [ = hm ( Wi{9) - Vb^0,i\ = (4.16) 

n— s>oo J J — ' n— s>oo \ ' / 



j=0 \ j=0 



We can now consider the dynamics. We can write out (4.6) in terms of 
our power series: 

oo ^ oo 



j=0 



-(c + £0(^)) J]£%(e,t) (4.17) 

j=0 



A formal treatment of the integrand could be considered, but without con- 
fidence of convergence, that seems unnecessary. We are only calculating uq 
and Ui, so we will ignore terms of 0(5^) or higher in the integrand: 

j=0 \j=0 U uo{uj,t)du) J Vj=o ^ ^ 



00 



+ £0(^)j^£%,(^,t) (4.18) 
j=0 

The above equation allows us to solve for the first two terms of the pertur- 
bation expansion. 

The equation for the first term uo{0,t) from the expansion is: 



(B^uoye,t) 



j uq{uj, t) auj 
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We can explicitly solve for the time-dependent total density by observing: 

d 



dt 

Solving the above gives: 



u{uj ^t) duj = 1 — / u{uj,t)duj 



u{u, t)du = 1 + Aexp[—t] 



(4.20) 



(4.21) 



where A = J u*{u)du. Substituting in to (4.19) gives: 



1 + Aexp[-t] 



(4.22) 



The right hand side of the above clearly depends upon the denominator 
1 + y4exp[— t] continuously in almost any operator topology. Since we are 
primarily concerned with the asymptotic dynamics, it is sufficient to show 
that the asymptotic equation converges: 



(4.23) 



From the results in Theorem [TJ we know that the above equation converges 
to a multiple of the principal eigenfunction wq. Equation (4.21) implies that 
the total density is equal to ^ and UQ{9,t) — )■ Wq{9) in norm as was required. 

The second term of the expansion {e^) is slightly more complicated: 



Juo{uj,t)duj [Juo{uj,t)dujy 



Substituting in known quantities gives: 

B-kUi j Ui 



Ui 



l + v4exp[-t] (1 + Aexp[-t])2 
We can now solve for the integral of Ui. 



B^uo){e,t)-u,{e,t)-mMo,t) 

(4.24) 
(4.25) 



B -k Uq — Ui — (f)Uo 



d_ 
dt 



Ui 



J Bi^Ui 

l + Aexp[-t] 1 + A exp [-t] 

- / Ml - 



Ml 



(4.26) 
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That gives the exphcit solution: 



ui{u, t) du = — exp 



t 

[-t] J exp[s] J 



(p{u)uo{uj, t) doj ds 



(4.27) 



Weak convergence of Uq — )■ Wq is sufficient to show that j Ui ~ j 'P- 
Substituting the asymptotic forms into (4.24) similar to the first term gives: 



ill = B i< ui — ui + 
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asymptotically. Finally, to show convergence, we can write Ui as: 

ui{d,t) = wi{d) + e{d,t) 



where e(^,0) = —Wi{6). Putting that form into (4.28) gives: 

e = i3 7^ e — e 



(4.28) 



(4.29) 



(4.30) 



Inspection of (4.13) shows that Wi{6) is bounded and is therefore in L^. 
Since B -k u — u has the Fourier eigenpairs {{■jn, exp[— m6'])}rigN as an 
orthonormal basis, we can write the norm of e as: 



^exp[7„t]| 



Wiin) 



(4.31) 



neN 



where Wi{n) is the n'th Fourier component of Wi{6). We know that 7„ < 
for all n > Cj^and sup„>i7„ < 0. Combining those facts with the fact that 
Wi{0) = because J wi = gives: 



^exp[7„t]|wi(n)|^ < exp 



n>l 



sup 7„ ) t 

n>l 



E 

n>l 



Wiin) 



(4.32) 



That implies that ui{6,t) — t- Wi{6) asymptotically exponentially. 

We have now shown that the first two components of the perturbation ex- 



pansion of (4.6) converge to an equilibrium. Calculating further terms would 



not provide any additional insight as the estimate from equation (4.10) gives 



a worse estimate after the ffist order. Recall from the proof of Proposition 



^the 7„'s are real since B is symmetric. 
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|4] that the integral of the equihbrium w has to be equal to the eigenvalue 
/i. However, it is easy to observe that the argument showing that J Wi = 
holds for all Wj, which implies that /X]j=o^"'^i ~ c ^'^^ ^' However, 
^^^Qe-'/ij 7^ /io = ^ for all k > 2. We can show that by caclulating the 
second term in the perturbation expansion. 

Using equation (4.10), we can gather the terms: 

-B-kW2- 'kWi + -kWo = H0W2 + HlWi + H2W0 (4.33) 

Observe that B-kWi = Wi + ^, and filling in other known quantities gives the 
relation: 

-{B-kW2- W2) = ^B -kWi-'^ + fi2Wo 

c 

B-kW2-W2 = ^ h — 

B^W2-W2 = ^ + !^ (4.34) 

C ZTT 

It is obvious that /i2 = <^==^ J (pwi = 0. However, we can observe that: 

/^ = /0(B-i)-^ = (0,(B-i)-: 

jB-I)i^i,(B-I)(B-I)-V' 
\b-1)wi,wi) (4.35) 



by observing that (B — l)wi = -. We can use the Cauchy-Schwarz inequality 
to show: 

{Bwi,wi) — {wi,wi) < ||Bwi||2||wi||2 — ||'U-'i||2 < (4.36) 

The last inequality comes from the fact that ||B/||2 < ||B||||/||2 = ||/||2 and 
that equality holds if and only if / = C. 

An explicit example can give insight into the non-convergence discussed 
in the previous paragraph. Assume B = ^ and k = l+e(f) where < 1. It 
is easy to see that the principle solution to (4.8) is w = and /i = / Y+i^- 
Moreover, we have the exact expansion in e: 

oo oo 

w = j2^'^j = J2'\-<py (4.37) 

j=0 j=0 
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By our hypothesis \e(f)\ < 1, we know the above converges. It is easy to see 
that: 

oo 

j=0 

which is the exact solution for the first term of the asymptotic expansion. 
However, the second term gives an incorrect solution: 



(B-I)- (0^-/0^) =Eb(/'^^-^ 

= j 02(cu) doo - (l)\e) ^ (l)\e) (4.39) 



Further terms would show the same difficulty. 
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All simulations were run using the Numpy [H [TTl [TTl [T8| extension to Python. 
The equations of motion were integrated using a simple Euler method. The 
simulations were run for 100 time units with a time step of 0.01 time units. 
The circle was discretized using 2^^ equally spaced points from — tt to vr (the 
power of 2 was used to speed up the fast Fourier transform). All integrals were 
taken using the trapezoidal method included in Numpy. The convolution 
was performed by taking the real fast Fourier transform of the branching 
kernel B and the density u{6, t), multiplying, and taking the inverse real fast 
Fourier transform. The built-in Numpy fast convolution method was not 
used because that method pads the two convolved functions with extra zeros 
to prevent circular convolution, but the equations used here explicitly call 
for the circular convolution. The convolution was normalized by dividing B 
by the integral of the convolution of B with the constant function B-k 
Finally, the total branching rate was normalized by integrating the density 
at the previous time step, i.e.: 



Using the integral from the previous time step and not a more sophisticated 
prediction-correction methods is justified by the following inequality: 



Since simulations remain bounded, the bound above can be made uniform. 



2126, 1996. 



A Methods 




(A.l) 
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The equilibrium was calculated a priori by iterating the equilibrium op- 
erator A'. Explicitly, a sequence of functions was generated by: 

Vn+m = r ]n^.A '^n{e) (A.3) 
j Vn[9) dO 

where the discretization and convolution were performed exactly as above 
and fo = ^ The theoretical justification for using this method is outlined 
in SectionlB} For three combinations of B and /t, ||fio4 — ^2x104111 was less 
than numerical precision. For the combination of Bi and K2-, fio^ was not 
sufficiently converged, so Vi^a was used. That decision was based on the 
condition that ||vio6 ~ 'y2xio'5||i was less than numerical precision. That level 
of precision was used to ensure that the convergence could be seen even 
when \\u{9,t) — Vn\\i < 10^^ — 10"^*^. Even 10^ iterations of A' only took 
several minutes on a standard Linux desktop system concurrently running 
other programs, a number that could be reduced with further optimization. 

At each time-step, the L} distance between the state of the system u{6, t) 
and the equilibrium was calculated. That quantity is plotted as a function 
of time in Figure [2] 



B Calculating Equilibrium Distributions 

The results in this article justify the use of a naive eigenvalue calculation algo- 
rithm. Calculating eigenvalues of integral equations is a non-trivial problem. 
Investigation into open questions regarding the generality of orientation pat- 
terns across branching and capping patterns, such as in [20], may require 
calculating the equilibrium solution to equations like the ones analyzed here. 
Moreover, in the previous section, equilibrium distributions were calculated 
a priori to show that simulations converged. The method below has proven 
to be very efficient for the work in this article. 

We will consider calculating the leading eigenvalue of the equilibrium 
operator for zeroth-order branching. A'. As A' is self-adjoint and compact, 
we can represent its range as the sum of eigenfunctions. We can explicitly 
calculate the n-th iterate of A' in terms of its (orthonormal) eigenfunctions: 

(a')"^; = (A')"(5^c,n,) = $^/i,"c,n,- (B.l) 

3 j 
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where Cj = {VyUj)^. We know that /ig is equal to the spectral radius from the 
Krein-Rutman theorem as in the proof of Theorem [T| The proof also implies 
the eigenvalue is simple. Finally, since A' is a compact operator, we know 
that there must be a spectral gap, i.e. fiQ — > c > for some c and all 

All that remains necessary to show that the above iteration converges to 
the positive equilibrium is to show that cj = (f,Mo)ft 7^ 0. If f equals the 
constant function, that condition is fulfilled. However, a stronger result is 
possible. By a result in [6], we know that Uq{6) > 0. The continuity of uq 
gives that inf Uq > 0. Thus, we have the inequality: 

{v,uo)^ > ml^{uo{9)K{9)) j v{e)de > (B.2) 
which implies that the iterative procedure will converge. 
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